1 Import data

languages <- read_csv("voicing-effect/stimuli/languages.csv")
## Parsed with column specification:
## cols(
##   speaker = col_character(),
##   language = col_character()
## )
words <- read_csv("voicing-effect/stimuli/nonce.csv")
## Parsed with column specification:
## cols(
##   item = col_integer(),
##   word = col_character(),
##   ipa = col_character(),
##   c1 = col_character(),
##   c1phonation = col_character(),
##   vowel = col_character(),
##   anteropost = col_character(),
##   height = col_character(),
##   c2 = col_character(),
##   c2phonation = col_character(),
##   c2place = col_character(),
##   language = col_character()
## )
columns <- c(
    "speaker",
    "seconds",
    "rec.date",
    "prompt",
    "label",
    "TT.displacement.sm",
    "TT.velocity",
    "TT.velocity.abs",
    "TD.displacement.sm",
    "TD.velocity",
    "TD.velocity.abs"
)

aaa.files <- list.files(
    path = "./voicing-effect/results/ultrasound",
    pattern = "*-tongue-cart.tsv",
    full.names = TRUE
)

tongues <- read_aaa(
    aaa.files,
    columns, 
    na.rm = TRUE
) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor) %>%
    group_by(speaker) %>%
    mutate(
        X.re = rescale(X),
        Y.re = rescale(Y)
    ) %>%
    ungroup() %>%
    mutate(
        vowel.ord = ordered(vowel, levels = c("a", "o", "u")),
        c2place.ord = ordered(c2place, levels = c("coronal", "velar")),
        c2phonation.ord = ordered(c2phonation, levels = c("voiceless", "voiced"))
    )
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character(),
##   X_1 = col_character(),
##   Y_1 = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Joining, by = "speaker"
## Joining, by = c("word", "language")
max <- tongues %>%
    filter(label %in% c("max_TT", "max_TD"), vowel != "u") %>%
    arrange(rec.date, fan.line) %>%
    create_event_start("rec.date")

max_it_12 <- max %>%
    filter(speaker %in% c("it01", "it02"))

max_pl_23 <- max %>%
    filter(speaker %in% c("pl02", "pl03"))

max_pl_34 <- max %>%
    filter(speaker %in% c("pl03", "pl04"))

aaa_files_clos <- list.files(
    path = "./voicing-effect/results/ultrasound",
    pattern = "*-tongue-clos-cart.tsv",
    full.names = TRUE
)

tongues_clos <- read_aaa(
    aaa_files_clos,
    columns, 
    na.rm = TRUE
) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor) %>%
    group_by(speaker) %>%
    mutate(
        X.re = rescale(X),
        Y.re = rescale(Y)
    ) %>%
    ungroup() %>%
    mutate(
        vowel.ord = ordered(vowel, levels = c("a", "o", "u")),
        c2place.ord = ordered(c2place, levels = c("coronal", "velar")),
        c2phonation.ord = ordered(c2phonation, levels = c("voiceless", "voiced"))
    ) %>%
    arrange(rec.date, fan.line) %>%
    create_event_start("rec.date")
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character(),
##   X_2 = col_character(),
##   Y_2 = col_character(),
##   X_3 = col_character(),
##   Y_3 = col_character(),
##   X_4 = col_character(),
##   Y_4 = col_character(),
##   X_5 = col_character(),
##   Y_5 = col_character(),
##   X_6 = col_character(),
##   Y_6 = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character(),
##   X_3 = col_character(),
##   Y_3 = col_character(),
##   X_4 = col_character(),
##   Y_4 = col_character(),
##   X_5 = col_character(),
##   Y_5 = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Joining, by = "speaker"
## Joining, by = c("word", "language")

2 Exploration

max %>%
    filter(c2place == "coronal", language == "italian") %>%
    plot_tongue(geom = "point", alpha = 0.5) +
    aes(colour = c2phonation) +
    facet_grid(vowel ~ speaker)

max %>%
    filter(c2place == "velar", language == "italian") %>%
    plot_tongue(geom = "point", alpha = 0.5) +
    aes(colour = c2phonation) +
    facet_grid(vowel ~ speaker)

max %>%
    filter(c2place == "coronal", language == "polish") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ speaker) +
    coord_fixed(1.5)

max %>%
    filter(c2place == "coronal", language == "polish") %>%
    ggplot(aes(X.re, Y.re, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(speaker ~ vowel) +
    coord_fixed(0.7)

max %>%
    filter(c2place == "velar", language == "polish") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ speaker) +
    coord_fixed(0.7)

max %>%
    filter(c2place == "velar", language == "polish") %>%
    ggplot(aes(X.re, Y.re, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(speaker ~ vowel) +
    coord_fixed(0.7)

3 Testing

Comparing fixed effects models works best with ML, otherwise you can use fREML with AIC for fixed effects, but if there is an AR1 model cannot use AIC, so need ML.

3.1 Polish

pl.m1 <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML"
)

rho <- start_value_rho(pl.m1)

pl.m1.ar <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML",
    rho = rho,
    AR.start = max_pl_23$start.event
)

summary(pl.m1)

pl.m1.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML"
)

compareML(pl.m1, pl.m1.null)
plot_smooth(
    pl.m1,
    view = "X.re",
    plot_all= "c2phonation.ord",
    cond = list("c2place.ord" = "coronal"),
    rug = FALSE
)
plot_diff(
    pl.m1,
    view = "X.re",
    comp = list(c2phonation.ord = c("voiceless", "voiced")),
    cond = list("c2place.ord" = "coronal")
)
plot_gamsd(
    pl.m1,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
plot_gamsd(
    pl.m1.ar,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)

3.1.1 PL02

pl02_max <- filter(max, speaker == "pl02", X > -20)

pl02.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl02.m1)

pl02.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "ML",
    rho = rho,
    AR.start = pl02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl02.m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.23973    0.59721   5.425 6.45e-08 ***
## X           -0.09449    0.03601  -2.624  0.00876 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df       F  p-value    
## s(X)                         7.269   7.773  66.456  < 2e-16 ***
## s(X):c2phonation.ordvoiced   8.003   8.684   6.916 5.04e-07 ***
## s(X):c2place.ordvelar        8.877   8.983 309.316  < 2e-16 ***
## s(X):vowel.ordo              7.803   8.573  21.311  < 2e-16 ***
## s(X,rec.date)              188.143 300.000   9.608  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 342/343
## R-sq.(adj) =  0.949   Deviance explained = 95.4%
## -ML =   5542  Scale est. = 4.4805    n = 2391
pl02.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl02.m1, pl02.m1.null)
## pl02.m1: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl02.m1.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##          Model    Score Edf  Chisq    Df   p.value Sig.
## 1 pl02.m1.null 5557.001  11                            
## 2      pl02.m1 5541.991  13 15.011 2.000 3.027e-07  ***
## 
## AIC difference: -35.85, model pl02.m1 has lower AIC.
plot_gamsd(
    pl02.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -19.967800 to 59.219800. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:14.

3.1.2 PL03

pl03_max <- filter(max, speaker == "pl03")

pl03.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl03.m1)

pl03.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML",
    rho = rho,
    AR.start = pl03_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl03.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML",
    rho = rho,
    AR.start = pl03_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl03.m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.15251    0.84915  -6.068 1.63e-09 ***
## X            0.42071    0.05257   8.004 2.35e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         7.433   7.793 52.044  < 2e-16 ***
## s(X):c2phonation.ordvoiced   6.740   7.654  4.817 1.54e-05 ***
## s(X):c2place.ordvelar        8.755   8.938 33.136  < 2e-16 ***
## s(X):vowel.ordo              8.770   8.941 20.139  < 2e-16 ***
## s(X,rec.date)              209.072 235.000 26.599  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 277/278
## R-sq.(adj) =  0.969   Deviance explained = 97.3%
## fREML = 3995.5  Scale est. = 2.8743    n = 1793
pl03.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl03.m1.ar, pl03.m1.ar.null)
## pl03.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl03.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl03.m1.ar.null 3335.866  11                         
## 2      pl03.m1.ar 3330.848  13 5.018 2.000   0.007  **
## Warning in compareML(pl03.m1.ar, pl03.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.611727, rho2 = 0.611727).
plot_gamsd(
    pl03.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): velar. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

3.1.3 PL04

filter(max, speaker == "pl04") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ c2place)

filter(max, speaker == "pl04", X > -30) %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ c2place)

pl04_max <- filter(max, speaker == "pl04", X > -30)

pl04.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl04.m1)

pl04.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "ML",
    rho = rho,
    AR.start = pl04_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl04.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "ML",
    rho = rho,
    AR.start = pl04_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl04.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.28296    0.38386   8.553   <2e-16 ***
## X           -0.03567    0.03469  -1.028    0.304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         7.538   7.847 31.856  < 2e-16 ***
## s(X):c2phonation.ordvoiced   1.006   1.009  0.296    0.589    
## s(X):c2place.ordvelar        8.539   8.876 26.159  < 2e-16 ***
## s(X):vowel.ordo              7.975   8.620 10.694 7.37e-15 ***
## s(X,rec.date)              197.731 230.000 25.617  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 272/273
## R-sq.(adj) =  0.975   Deviance explained =   98%
## -ML = 1757.9  Scale est. = 0.72325   n = 1228
pl04.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
# compareML(pl04.m1, pl04.m1.null)
compareML(pl04.m1.ar, pl04.m1.ar.null)
## pl04.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl04.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl04.m1.ar.null 1758.008  11                         
## 2      pl04.m1.ar 1757.861  13 0.148 2.000   0.863
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.518140, rho2 = 0.518140).
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): Only small difference in ML...
plot_gamsd(
    pl04.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -29.959700 to 32.556800. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 26/05/2017 17:58:11.

ggsave("presentations/2017-ultrafest/pl04.pdf", width = 6, height = 4)
acf_plot(resid(pl04.m1.ar), split_by=list(pl04_max$rec.date))

acf_resid(pl04.m1.ar, split_pred = "AR.start")

pl04_a <- filter(max, speaker == "pl04", vowel == "a")

pl04.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "fREML"
)

rho <- start_value_rho(pl04.m1)

pl04.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "ML",
    rho = rho,
    AR.start = pl04_a$start.event
)

pl04.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "ML",
    rho = rho,
    AR.start = pl04_a$start.event
)

summary(pl04.m1.ar)
compareML(pl04.m1.ar, pl04.m1.ar.null)

If only /a/ is used, no significance at all. With /o/ and X > -30, no significance, but with all X some significant difference at the very left end of the tongue.

plot_gamsd(
    pl04.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)

3.1.4 PL03-PL04

pl34.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "fREML"
)

rho <- start_value_rho(pl34.m1)

pl34.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)

pl34.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)

summary(pl34.m1.ar)
compareML(pl34.m1.ar, pl34.m1.ar.null)
plot_gamsd(
    pl34.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar", speaker = "pl04")
)

The models with X, Y are not reliable because the tongue size differs a lot between speakers.

pl34.m1.re <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "fREML"
)

rho <- start_value_rho(pl34.m1.re)

pl34.m1.ar.re <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)

pl34.m1.ar.re.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)

summary(pl34.m1.ar.re)
compareML(pl34.m1.ar.re, pl34.m1.ar.re.null)
plot_gamsd(
    pl34.m1.ar.re,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
plot_gamsd(
    pl34.m1.ar.re,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar")
)

3.2 Italian

it01 and it02 have TRA at maximum displacement.

it.m1 <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "fREML"
)

rho <- start_value_rho(it.m1)

it.m1.ar <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML",
    rho = rho,
    AR.start = max_it_12$start.event
)

summary(it.m1.ar)

it.m1.ar.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
#        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML",
    rho = rho,
    AR.start = max_it_12$start.event
)

it.m1.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML"
)

compareML(it.m1.ar, it.m1.ar.null)
plot_gamsd(
    it.m1.ar,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)

3.2.1 IT01

it01_max <- filter(max, speaker == "it01")

it01.m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it01.m1)

it01.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "ML",
    rho = rho,
    AR.start = it01_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it01.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.73141    0.30682  -15.42   <2e-16 ***
## X            0.76532    0.02666   28.70   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df       F  p-value    
## s(X)                        7.764   7.944 181.554  < 2e-16 ***
## s(X):c2phonation.ordvoiced  5.037   6.136   9.480 2.06e-10 ***
## s(X):c2place.ordvelar       8.822   8.979 183.531  < 2e-16 ***
## s(X):vowel.ordo             6.862   7.808  11.566 1.19e-15 ***
## s(X,rec.date)              92.680 225.000   1.976  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 267/268
## R-sq.(adj) =  0.979   Deviance explained =   98%
## -ML = 3488.2  Scale est. = 4.0701    n = 1932
it01.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "ML",
    rho = rho,
    AR.start = it01_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it01.m1.ar, it01.m1.ar.null)
## it01.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it01.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf  Chisq    Df   p.value Sig.
## 1 it01.m1.ar.null 3509.547  11                            
## 2      it01.m1.ar 3488.231  13 21.316 2.000 5.527e-10  ***
## Warning in compareML(it01.m1.ar, it01.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.745137, rho2 = 0.745137).
plot_gamsd(
    it01.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -40.731200 to 54.425500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:10:52.

ggsave("presentations/2017-ultrafest/it01.pdf", width = 6, height = 4)

3.2.2 IT02

it02_max <- filter(max, speaker == "it02")

it02.m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it02.m1)

it02.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "ML",
    rho = rho,
    AR.start = it02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it02.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.67360    0.52437  -3.192  0.00145 ** 
## X            0.65159    0.03157  20.640  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df      F  p-value    
## s(X)                        6.386   7.189 57.403  < 2e-16 ***
## s(X):c2phonation.ordvoiced  6.589   7.613  4.927 9.93e-06 ***
## s(X):c2place.ordvelar       8.630   8.910 52.382  < 2e-16 ***
## s(X):vowel.ordo             5.937   7.020  9.062 6.16e-11 ***
## s(X,rec.date)              72.997 185.000  1.401  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 227/228
## R-sq.(adj) =  0.916   Deviance explained = 92.2%
## -ML = 2803.6  Scale est. = 9.778     n = 1233
it02.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "ML",
    rho = rho,
    AR.start = it02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it02.m1.ar, it02.m1.ar.null)
## it02.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it02.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df   p.value Sig.
## 1 it02.m1.ar.null 2813.488  11                           
## 2      it02.m1.ar 2803.592  13 9.897 2.000 5.034e-05  ***
## Warning in compareML(it02.m1.ar, it02.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.735050, rho2 = 0.735050).
plot_gamsd(
    it02.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -44.949000 to 62.666500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:45:14.

4 Italian tongue at closure

tongues_clos %>%
    filter(c2place == "coronal", language == "italian") %>%
    plot_tongue(geom = "point", alpha = 0.5) +
    aes(colour = c2phonation) +
    facet_grid(vowel ~ speaker)

4.1 IT01

it01_clos <- tongues_clos %>%
    filter(speaker == "it01", vowel != "u")
it01_clos_m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_clos,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it01_clos_m1)

it01_clos_m1_ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_clos,
    method = "ML",
    rho = rho,
    AR.start = it01_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it01_clos_m1_ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61066    0.39198   1.558    0.119    
## X            0.84888    0.02509  33.831   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df       F  p-value    
## s(X)                        7.549   7.867 198.722  < 2e-16 ***
## s(X):c2phonation.ordvoiced  4.683   5.594  13.824 3.56e-14 ***
## s(X):c2place.ordvelar       8.641   8.924 174.974  < 2e-16 ***
## s(X):vowel.ordo             6.049   6.977  14.770  < 2e-16 ***
## s(X,rec.date)              96.251 230.000   1.789  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 272/273
## R-sq.(adj) =  0.981   Deviance explained = 98.2%
## -ML = 2983.9  Scale est. = 3.7241    n = 1648
it01_clos_m1_ar_null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_clos,
    method = "ML",
    rho = rho,
    AR.start = it01_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it01_clos_m1_ar, it01_clos_m1_ar_null)
## it01_clos_m1_ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it01_clos_m1_ar_null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                  Model    Score Edf  Chisq    Df   p.value Sig.
## 1 it01_clos_m1_ar_null 3014.072  11                            
## 2      it01_clos_m1_ar 2983.854  13 30.218 2.000 7.523e-14  ***
## Warning in compareML(it01_clos_m1_ar, it01_clos_m1_ar_null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.721734, rho2 =
## 0.721734).
plot_gamsd(
    it01_clos_m1_ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -39.125800 to 54.502000. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:10:52.

ggsave("presentations/2017-ultrafest/it01-clos.pdf", width = 6, height = 4)

4.2 IT02

it02_clos <- tongues_clos %>%
    filter(speaker == "it02", vowel != "u")
it02_clos_m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_clos,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it02_clos_m1)

it02_clos_m1_ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_clos,
    method = "ML",
    rho = rho,
    AR.start = it02_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it02_clos_m1_ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.56479    0.48605  -5.277 1.57e-07 ***
## X            0.64094    0.04079  15.714  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         7.047   7.612 40.215  < 2e-16 ***
## s(X):c2phonation.ordvoiced   7.191   8.131  5.061 2.41e-06 ***
## s(X):c2place.ordvelar        8.088   8.702 29.725  < 2e-16 ***
## s(X):vowel.ordo              8.257   8.776 10.910 2.49e-15 ***
## s(X,rec.date)              115.667 195.000  2.711  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 237/238
## R-sq.(adj) =   0.92   Deviance explained = 92.8%
## -ML = 2848.9  Scale est. = 6.7751    n = 1319
it02_clos_m1_ar_null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_clos,
    method = "ML",
    rho = rho,
    AR.start = it02_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it02_clos_m1_ar, it02_clos_m1_ar_null)
## it02_clos_m1_ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it02_clos_m1_ar_null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                  Model    Score Edf  Chisq    Df   p.value Sig.
## 1 it02_clos_m1_ar_null 2861.022  11                            
## 2      it02_clos_m1_ar 2848.944  13 12.078 2.000 5.680e-06  ***
## Warning in compareML(it02_clos_m1_ar, it02_clos_m1_ar_null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.718876, rho2 =
## 0.718876).
plot_gamsd(
    it02_clos_m1_ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): o. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:44:52.

Both it01 and it02 show TRA at closure. It is difficult to quantify the difference compared to TRA at maximum displacement. Need a model to check that.

5 Polish tongue at closure

5.1 PL02

pl02_clos <- tongues_clos %>%
    filter(speaker == "pl02", vowel != "u", X > -21)
pl02_clos_m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_clos,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl02_clos_m1)

pl02_clos_m1_ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_clos,
    method = "ML",
    rho = rho,
    AR.start = pl02_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl02_clos_m1_ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.04534    0.36324  13.890   <2e-16 ***
## X           -0.03225    0.03627  -0.889    0.374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df       F  p-value    
## s(X)                         6.593   7.454  16.385  < 2e-16 ***
## s(X):c2phonation.ordvoiced   1.003   1.004   0.301    0.584    
## s(X):c2place.ordvelar        8.724   8.942 122.679  < 2e-16 ***
## s(X):vowel.ordo              4.633   5.716   4.795 9.78e-05 ***
## s(X,rec.date)              168.762 315.000   5.365  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 357/358
## R-sq.(adj) =  0.908   Deviance explained = 91.7%
## -ML = 3963.2  Scale est. = 4.0457    n = 1897
pl02_clos_m1_ar_null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_clos,
    method = "ML",
    rho = rho,
    AR.start = pl02_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl02_clos_m1_ar, pl02_clos_m1_ar_null)
## pl02_clos_m1_ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl02_clos_m1_ar_null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                  Model    Score Edf Chisq    Df p.value Sig.
## 1 pl02_clos_m1_ar_null 3963.339  11                         
## 2      pl02_clos_m1_ar 3963.189  13 0.150 2.000   0.860
## Warning in compareML(pl02_clos_m1_ar, pl02_clos_m1_ar_null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.514127, rho2 =
## 0.514127).
## Warning in compareML(pl02_clos_m1_ar, pl02_clos_m1_ar_null): Only small difference in ML...
plot_gamsd(
    pl02_clos_m1_ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -20.995600 to 42.515500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:21:39.

5.2 PL04

pl04_clos <- tongues_clos %>%
    filter(speaker == "pl04", vowel != "u")
pl04_clos_m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_clos,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl04_clos_m1)

pl04_clos_m1_ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_clos,
    method = "ML",
    rho = rho,
    AR.start = pl04_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl04_clos_m1_ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.53763    0.11242 -13.677   <2e-16 ***
## X           -0.01176    0.01746  -0.674    0.501    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df       F  p-value    
## s(X)                         7.817   7.951 234.521  < 2e-16 ***
## s(X):c2phonation.ordvoiced   1.004   1.005   0.809    0.369    
## s(X):c2place.ordvelar        5.960   7.099   5.674 1.66e-06 ***
## s(X):vowel.ordo              7.493   8.406  16.701  < 2e-16 ***
## s(X,rec.date)              155.677 234.000   5.321  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 277/278
## R-sq.(adj) =  0.984   Deviance explained = 98.6%
## -ML = 1886.4  Scale est. = 0.83789   n = 1387
pl04_clos_m1_ar_null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_clos,
    method = "ML",
    rho = rho,
    AR.start = pl04_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl04_clos_m1_ar, pl04_clos_m1_ar_null)
## pl04_clos_m1_ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl04_clos_m1_ar_null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                  Model    Score Edf Chisq    Df p.value Sig.
## 1 pl04_clos_m1_ar_null 1886.770  11                         
## 2      pl04_clos_m1_ar 1886.368  13 0.401 2.000   0.669
## Warning in compareML(pl04_clos_m1_ar, pl04_clos_m1_ar_null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.494196, rho2 =
## 0.494196).
## Warning in compareML(pl04_clos_m1_ar, pl04_clos_m1_ar_null): Only small difference in ML...
plot_gamsd(
    pl04_clos_m1_ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -39.613600 to 30.560600. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 26/05/2017 17:59:34.

ggsave("presentations/2017-ultrafest/pl04-clos.pdf", width = 6, height = 4)

6 Comparison of Italian voiced at closure and maximum displacement

6.1 IT01

it01_voiced <- rbind(it01_max, it01_clos) %>%
    filter(c2phonation == "voiced") %>%
    mutate(
        position = ifelse(label %in% c("max_TT", "max_TD"), "maximum", "closure"),
        position_ord = ordered(position, levels = c("maximum", "closure"))
    ) %>%
    unite(item_no, seconds:rec.date) %>%
    filter(vowel != "u") %>%
    mutate_if(is.character, as.factor)
it02_voiced_m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = position_ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, item_no, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_voiced,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it02_voiced_m1)

it02_voiced_m1_ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = position_ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, item_no, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_voiced,
    method = "ML",
    rho = rho,
    AR.start = it01_voiced$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it02_voiced_m1_ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = position_ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, item_no, bs = "fs", xt = "cr", m = 1, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.06151    0.32749  -6.295 3.92e-10 ***
## X            0.96067    0.02656  36.166  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                             edf  Ref.df       F  p-value    
## s(X)                      7.594   7.881 313.331  < 2e-16 ***
## s(X):position_ordclosure  3.981   4.904  11.740 6.68e-11 ***
## s(X):c2place.ordvelar     8.186   8.585 219.239  < 2e-16 ***
## s(X):vowel.ordo           6.676   7.561  16.233  < 2e-16 ***
## s(X,item_no)             79.093 225.000   1.544  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 267/268
## R-sq.(adj) =   0.98   Deviance explained = 98.1%
## -ML = 3204.3  Scale est. = 3.8788    n = 1776
it02_voiced_m1_ar_null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, item_no, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_voiced,
    method = "ML",
    rho = rho,
    AR.start = it01_voiced$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it02_voiced_m1_ar, it02_voiced_m1_ar_null)
## it02_voiced_m1_ar: Y ~ X + s(X, bs = "cr") + s(X, by = position_ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, item_no, bs = "fs", xt = "cr", m = 1, k = 5)
## 
## it02_voiced_m1_ar_null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, item_no, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                    Model    Score Edf  Chisq    Df   p.value Sig.
## 1 it02_voiced_m1_ar_null 3224.556  11                            
## 2      it02_voiced_m1_ar 3204.308  13 20.248 2.000 1.608e-09  ***
## Warning in compareML(it02_voiced_m1_ar, it02_voiced_m1_ar_null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.727873, rho2 =
## 0.727873).
plot_gamsd(
    it02_voiced_m1_ar,
    view = "X",
    comparison = list(position_ord = c("maximum", "closure")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -39.125800 to 53.346500. 
##  * position_ord : factor; set to the value(s): closure, maximum. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * item_no : factor; set to the value(s): 1.0035_29/11/2016 15:16:05.

ggsave("presentations/2017-ultrafest/it01-voiced.pdf", width = 6, height = 4)